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13  ABSTRACT  ^  v  .T.) Forecast  errors  exhibit  the  characteristics  of  approximations 

in  simulating  dynamical  and  physical  processes  in  models.  The  models  are  very  com¬ 
plex  and  hence  it  is  not  always  possible  to  identify  the  approximations  responsible 
for  any  particular  error  pattern  in  forecasts.  A  comparison  between  the  models' 
forecast  performances  can  be  valuable  in  isolating  the  causes  of  error  patterns. 

Here  a  comparison  of  forecast  errors  in  the  AFGL  and  COLA  models  is  made  with  the 
intent  of  identifying  the  causes  of  forecast  errors.  The  two  models  are  based  on 
identical  approximations  in  simulating  the  dynamical  processes  and  only  minor  differ 
ences  in  parameterlzatlons  of  the  physical  processes.  Nine  ten-day  forecasts  are 
made  to  study  the  error  characteristics  in  the  two  models.  The  errors  in  the  500  mb 
geopotential  height  are  negative  in  tropics  and  positive  in  extratropics.  The  tem¬ 
peratures  at  850  mb  are  colder  than  observed  in  tropics  and  warmer  than  observed  in 
extratropics.  At  150  mb  the  temperatures  are  warmer  than  observed  in  tropics  and 
colder  than  observed  in  extratropics.  These  qualitative  error  characteristics  are 
not  only  common  to  these  two  models,  but  also  they  are  common  in  the  NMC,  GFDL  and 
ECMWF  forecast  models.  The  difference  in  the  error  structure  between  the  two  models 
is  the  magnitude  of  the  error  in  tropics.  The  tropical  error  in  the  AFGL  model  is 
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larger  chan  that  In  the  COLA  model.  Another  difference  is  In  the  830  mb  relative 
humidity  field.  In  the  AFGL  model, relative  humidity  errors  are  negative  largely 
over  the  ocean  and  positive  over  land  with  minor  exceptions.  This  error  structure 
differs  from  that  of  the  COLA  model  which  consists  of  mostly  positive  errors  every¬ 
where  with  some  small  regions  of  negative  errors.  The  major  differences  in  the 
physical  parameterlzations  between  the  two  models  is  in  the  radiation  interaction 
with  deep  convective  clouds,  the  manner  in  which  the  sea  surface  temperature  (SST) 
is  prescribed  and  the  vertical  transport  of  heat  and  moisture  by  shallow  convection. 
The  magnitude  of  tropical  errors  in  the  geopotential  height  is  500  mb  and  temperature 
at  850  ^  may  be  because  the  AFGL  model  does  not  include  deep  convective  cloud- 
.adlatlon  interactions.  The  850  mb  relative  humidity  errors  over  oceans  are  probably 
due  to  Che  manner  in  which  the  SST  is  prescribed  and  the  lack  of  proper  vertical 
transport  of  moisture  by  the  shallow  convection  parameterization. 
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Preface 


The  purpose  of  the  research  carried  out  under  the  contract  No.  F19628-88-K-0015 
was  to  determine  the  impact  of  orographic  effects  on  the  short  and  medium-range  forecasts 
with  the  AFGL-GSM.  The  orographic  effects  include  large  scale  blocking  effects  of 
massive  mountain  ranges  over  the  globe  and  the  subgrid  scale  effects  of  mountains  of  width 
of  about  50  km  on  the  large  scale  flow.  The  large  scale  blocking  effects  are  represented  by 
the  silhouette  orography.  The  subgrid  scale  effects  are  the  effects  of  orographic  gravity 
wave  drag.  The  gravity  waves  aid  in  transferring  momentum  between  the  earth’s  surface 
and  the  atmosphere.  This  exchange  could  be  in  all  levels  of  the  troposphere  and  in  the 
stratosphere.  The  impact  of  the  orographic  effects  and  of  horizontal  resolution  on  the  short 
and  medium  range  forecasts  were  studied  on  an  earlier  version  of  the  AFGL  model.  The 
results  were  reported  in  Zhou  (1990).  The  parameterization  of  gravity  wave  drag  in  this 
study  was  based  on  a  linear  theory.  The  impact  of  gravity  wave  drag  parameterization 
based  on  nonlinear  theory  was  reported  in  Kirtman  et  al.  (1991).  Now  the  silhouette 
orography  and  the  gravity  wave  drag  parameterization  based  on  nonlinear  theory  are 
implemented  in  the  current  version  of  the  AFGL  model.  The  gravity  wave  drag 
parameterization  was  also  implemented  in  the  COLA  (Center  for  Ocean-Land-Atmosphere 
Interactions)  model.  The  purpose  of  this  report  is  to  present  a  comparison  of 

y  ■*%  •  •  ^ 

— medij^ni-jange  forecast  performances  of  the  AFGL  model  with  that  of  the  COLA  model 
with  th.e  of  identifying  the  deficiencies  in  simulating  dynamical  and  physical 
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1.  INTRODUCTION 

Accuracy  of  numerical  weather  prediction  is  limited  by  the  approximations  in  the 
predictive  system  and  the  errors  in  determining  the  initial  state  of  the  atmosphere.  In  a 
perfect  predictive  system  the  accuracy  of  the  forecast  would  deteriorate  with  time  due  to 
nonlinear  interactions  (Lorenz,  1969)  and  hydrodynamical  instabilities  (Leith,  1971)  so  long 
as  there  is  a  non-zero  error  in  the  initial  state.  The  rate  at  which  the  forecast  deteriorates 
depends  on  the  growth  rate  of  instabilities,  the  nature  of  nonlinear  interactions  and  the 
amplitude  and  structure  of  the  initial  error.  The  total  forecast  error  is  the  difference 
between  the  forecast  and  the  observation  and  is  commonly  measured  by  the  difference 
between  the  forecast  and  the  analysis  at  the  verification  time.  The  total  forecast  error  is 
therefore  due  to  both  the  imperfections  in  the  predictive  system  and  the  error  in  the  initial 
state.  An  ensemble  average  of  the  total  forecast  errors  based  on  a  large  number  of  forecasts 
made  with  synoptically  independent  initial  conditions  is  referred  to  as  the  systematic  error. 
For  any  particular  day’s  forecast  the  difference  between  the  mean  square  of  the  total  error 
and  the  mean  square  of  the  systematic  error  is  referred  to  as  the  mean  square  of  the 
transient  error.  The  systematic  error  is  usually  considered  to  be  due  to  imperfections  in 
the  model;  however,  the  systematic  errors  and  the  transient  errors  are  not  decoupled  in  any 
model. 

A  large  number  of  studies  have  been  made  to  identify  the  spatial  structure  of 
systematic  errors  (see,  eg.,  Fawcett,  1969;  Hollingsworth  et  al.  1980;  Wallace  and 
Woessner,  1981;  Arpe  and  Klinker,  1986;  among  others).  The  systematic  errors  in  models 
usually  have  coherent  patterns  in  space.  These  geographically  fixed  error  patterns  largely 
change  in  amplitude  with  the  increase  in  forecast  time.  The  coherent  structures  of  the 
systematic  errors  appear  in  zonally  symmetric  fields  as  well  as  in  the  zonally  asymmetric 
fields.  In  a  comprehensive  study  comparing  wintertime  systematic  errors  in  the  ECMWF 
(European  Centre  for  Medium  Range  Weather  Forecasts)  and  the  GFDL  (Geophysical 


1 


Fluid  Dynamics  Laboratory)  forecast  models,  Hollingsworth  et  al.  (1980)  estimated  that 
the  systematic  error  accounts  for  more  than  one  third  of  the  total  error  in  medium-range 
forecasts.  Hence,  a  significant  improvement  can  be  made  in  the  forecast  skill  by  reducing 
the  systematic  error.  Assuming  that  the  growth  rate  of  systematic  error  was  linear 
Miyakoda  et  al.  (1986)  suggested  that  the  forecast  can  be  improved  by  subtracting  the 
systematic  error  patterns  from  the  forecast  fields.  Schemm  and  Faller  (1986)  proposed  a 
scheme  to  statistically  correct  the  forecast  during  integration  from  empirical  relationships 
between  the  model  variables  and  the  systematic  errors.  These  procedures  are  only  stop  gap 
measures  before  we  identify  the  deficiencies  in  the  representations  of  dynamical  and 
physical  processes  responsible  for  the  systematic  errors  and  improve  the  representations  of 
these  processes. 

Bettge  (1983)  showed  that  the  systematic  errors  in  the  500  mb  geopotential  height 
forecasts  from  ECMWF  and  NMC  (National  Meteorological  Center)  models  have  similar 
patterns  in  that,  negative  error  patterns  are  observed  over  the  Rockies,  the  Alps  and  the 
Himalayan  mountain  ranges.  He  argued  that  the  negative  errors  are  due  to  the  model 
orography  because  it  is  computed  by  averaging  over  the  model  grid  which  underestimates 
the  height  of  high  mountains  He  considered  the  conservation  of  potential  vorticily  of  the 
flow  when  passing  over  the  model  orography.  As  the  vertical  extent  of  the  atmosphere  is 
decreased  over  the  top  of  mountain  an  anticyclonic  vorticity  is  generated  to  conserve  the 
potential  vorticity.  If  the  model  mountains  are  shallow  compared  to  that  of  the  earth,  the 
model  will  not  generate  as  much  anticyclonic  vorticity  and  that  will  result  in  the  negative 
errors.  Wallace  et  al  (1983)  enhanced  the  model  orography  by  adding  a  constant  multiple 
of  the  orographic  standard  deviation  to  the  mean  orography.  The  standard  deviation  of 
orography  was  computed  from  a  high  resolution  orography  data.  They  showed  that  the 
enhanced  orography  in  the  model  reduces  the  systematic  error  in  the  vicinity  of  the 
mountain  ranges. 
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Recently,  when  horizontal  and  vertical  resolutions  were  increased,  new  systematic 
errors  were  identified  which  were  absent  in  the  low  resolution  models.  The  errors  were 
positive  in  the  midlatitudo  zonally-averaged  zonal  wind.  Also,  there  were  negative  errors 
in  the  high-latitudes  zonally-averaged  temperature  especially  in  northern  hemisphere 
winter.  In  the  low  resolution  models  both  the  surface  drag  over  the  mountains  and  the 
horizontal  momentum  flux  convergence  into  the  middle  and  high  latitudes  in  the  upper  and 
lower  stratosphere  by  eddies  were  underestimated.  In  these  models,  the  simulated  northern 
hemispheric  westerly  jet  agreed  well  with  observations  but  the  southern  hemispheric  jet 
was  slightly  weak  When  the  horizontal  resolution  increased  the  effects  of  explicitly 
resolved  eddies  was  to  increase  the  horizontal  momentum  flux  convergence.  This  increased 
the  strength  of  the  westerly  jet  in  the  northern  hemisphere  but  the  southern  hemispheric 
westerly  jet  was  in  good  agreement  with  observations.  The  cold  temperature  bias  in  high 
latitudes  was  consistent  with  the  positive  westerly  bias  to  approximately  maintain  the 
thermal  wind  relation  One  way  of  reducing  the  positive  westerly  bias  is  to  increase  the 
dissipative  processes  of  zonal  momentum  in  the  northern  hemisphere  winter.  Lilly  (1972) 
proposed  that  the  dissipative  role  of  orographically  induced  gravity  waves  was  important 
enough  to  include  their  effects  explicitly  in  the  numerical  weather  prediction  models  and 
the  general  circulation  models.  A  large  number  of  studies  have  now  shown  that  the  explicit 
representations  of  the  effect  of  gravity  wave  drag  in  the  model  reduce  the  westerly  bias  in 
the  zonal  wind  and  the  cidd  bias  in  the  temperature  (see  e.g.,  Palmer  et  al.  (1989);  Helfand 
et  al.  (1987);  McFarlane  ( 1 987);  Iwasaki  et  al.  (1989);  Kirtman  et  al.  (1991)) 

Determining  the  spatial  structure  of  the  systematic  error  field  is  a  straightforward 
procedure;  however,  it  is  not  obvious  how  to  associate  the  error  field  with  the  physical  or 
dynamical  approximations  in  the  model.  A  comparison  between  the  results  of  different 
models  as  they  relate  to  the  real  atmosphere  can  provide  us  with  clues  to  isolate  the  causes 
of  errors  so  that  further  improvements  in  the  parameterizations  can  be  made  to  reduce  the 
forecast  errors. 
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The  purpose  of  this  study  is  to  compare  the  systematic  errors  of  the  AFGL  and 
COLA  forecast  models  with  the  intent  of  isolating  the  causes  of  errors  so  that  they  can  be 
reduced  by  improved  parameterization.  The  next  section  gives  a  brief  comparison  of  the 
two  models’  dynamical  and  physical  parameterizations.  The  third  section  explains  the 
experiments  and  analysis  procedure.  The  fourth  section  discusses  the  results  A  summary 
and  conclusions  are  presented  in  the  fifth  section. 

2  FORECAST  MODELS 

Both  the  AFGL  and  COLA  forecast  models  owe  their  numerical  frame  of  the 
dynamics  to  the  NMC  spectral  model  of  Sela  (1980).  Therefore,  the  simulation  of 
atmospheric  dynamical  processes  are  identical.  The  large  scale  as  well  as  subgrid  scale 
orographic  effects  are  treated  in  the  same  manner  in  both  models.  Both  models  include 
physical  parameterizations  of  planetary  boundary  layer,  shallow  convection,  deep  moist 
convection,  diurnal  variations  and  cloud-radiation  interactions  The  only  difference 
between  the  two  models  is  the  representations  of  these  physical  parameterizations  Here 
we  shall  describe  briefly  the  treatment  of  dynamics  and  differences  in  the  representations  of 
the  physical  processes. 

The  dynamics  of  the  model  is  based  on  the  primitive  equations  The  prognostic 
equations  prescribe  the  time  changes  of  vorticity,  divergence,  thermodynamic  energy, 
continuity  of  water  vapor  and  natural  log  of  surface  pressure  Vertical  a-velocity  and 
geopotential  height  are  computed  from  diagnostic  equations  The  variables  are  represented 
by  spherical  harmonics  with  rhomboidal  truncation  at  wave  number  dO  The  horizontal 
transform  grid  is  equally  spaced  in  the  east-west  direction  with  96  grid  points  along  the 
latitude  circle.  In  the  north-south  direction,  there  are  80  points  on  the  Gausian  grid.  The 
horizontal  resolution  is  therefore  approximately  2  25°  x  3.75°  lat-long  grid  Vertical 
coordinate  is  a  =  p/p*  where  p*  is  the  surface  pressure  There  are  19  levels  in  the  vertical 
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The  levels  are  spaced  in  vertical  as  1.00,  0.99,  0.973,  0  948,  0.893,  0.82,  0.735,  0.642,  0.546, 
0.45,  0.40,  0.35,  0  30,  0.25,  0.20,  0.15,  0.10,  0.05,  0.00  in  a.  The  semi-implicit  integration 
is  used  with  15  minute  time  interval.  The  horizontal  diffusion  in  the  model  is  scale 
selective  that  is,  where  F  is  a  prognostic  variable.  The  values  of  the  diffusion 

constant,  k,  and  the  manner  in  which  the  diffusion  term  is  used  for  the  prognostic  variables 
differs  in  the  two  models.  In  the  AFGL  model  diffusion  term  is  applied  to  all  modes  of  the 
divergence  but  for  the  vorticity,  temperature  and  specific  humidity  it  is  applied  only  to 
modes  in  the  upper  half  of  the  rhomboid.  The  magnitude  of  the  diffusion  constant  is  the 
same  in  both  cases,  which  is  k  =  6  x  lO*^  nr*  sec*.  In  the  COLA  model  k  is  computed  from 
k  =  where  a  is  the  radius  of  the  earth.  N  is  the  highest  wave  number 

resolved,  that  is,  N  =  3('  in  this  model,  t  is  the  dissipation  time  in  seconds.  For  the 
divergence  t  is  21  minutes  and  for  the  vorticity,  temperature  and  specific  humidity  r  is  28 
minutes.  Hence,  the  value  of  k  is  2.5  x  10>®  for  the  divergence  and  1.9  x  10''’  for  the 
vorticity,  temperature  and  specific  humidity. 

The  large  scale  orography  effects  are  simulated  by  the  silhouette  orography.  It  is 
computed  from  the  global  lO-minute  elevation  data  from  the  U.S.  Navy.  At  any  grid  point 
the  silhouette  orography  is  computed  as  an  average  of  maximum  peaks  in  profiles  of 
mountains  in  the  east-west  and  north-south  directions  over  the  grid  box  (see.  Appendix 
A).  Silhouette  orography  enhances  the  topography  in  a  similar  manner  as  in  the  envelope 
orography  but  yields  a  smoulher  surface  than  the  envelope  orography  especially  in  regions 
of  large  subgrid  scale  variance  in  a  narrow  zone,  along  the  southern  edge  of  Tibetan 
plateau,  for  example.  There  is  a  slight  difference  between  the  silhouette  orographies  in  the 
two  models  In  the  .XFCL  model  the  spectral  form  of  the  orography  was  smoothed  (see, 
.•\ppendix  a)  to  reduce  the  amplitude  of  the  Gibbs  oscillations  over  ocean.  In  the  COLA 
model  no  smoothing  was  applied. 

Depending  on  the  atmospheric  stability,  the  gravity  waves  induced  by  the  subgrid 
scale  orography  prop, agate  vertically  until  they  are  dissipated  in  the  critical  layer  or  as  a 
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consequence  of  convective  and  shear  instabilities.  These  waves  transport  horizontal 
momentum  upward  and  deposit  this  momentum  where  the  waves  are  ultimately  dissipated. 
The  momentum  is  lost  by  the  large  scale  where  the  waves  are  dissipated  and  it  is  brought 
down  to  the  earth  surface  where  it  is  deposited.  The  parameterization  of  these  effects 
consist  of  surface  wave  drag  and  the  vertical  distribution  of  the  wave  drag.  The  surface 
wave  drag  is  parameterized  by  following  the  procedure  suggested  by  Pierhumbert  (1987) 
and  the  vertical  variation  is  parameterized  according  to  Palmer  et  al.  (1986)  and  Helfand 
et  al.  (1987)  (see,  Appendix  B  for  methodology  and  Appendix  C  for  computer  code). 

The  atmosphere  exchanges  momentum,  heat  and  water  vapor  with  the  earth’s  surface 
in  the  planetary  boundary  layer.  This  exchange  depends  on  the  surface  cover.  Over  land, 
the  earth  may  be  covered  by  snow,  ice,  and  a  variety  of  different  plants,  sands  and  soils. 
Over  the  ocean,  the  surface  characteristics  are  determined  by  sea  ice  and  spatially  varying 
sea  surface  temperature.  In  the  AFGL  model,  the  planetary  boundary  over  land  consists  of 
three  parts;  soil  layer,  surface  layer  and  turbulent  mixed  layer.  The  soil  layer  exchanges 
heat  and  water  vapor  with  the  surface  layer.  The  surface  layer  exchanges  heat,  momentum 
and  water  vapor  with  the  planetary  boundary  layer.  In  the  soil  layer,  heat  is  transferred 
by  diffusion  and  water  by  diffusion  and  gravitational  transport  (see,  Mahrt  and  Pan,  1984). 
Surface  fluxes  of  heat,  momentum  and  water  vapor  are  represented  by  similarity  theory 
according  to  Louis  (1979).  The  turbulent  mixed  layer  grows  in  height  due  to  the  surface 
heating  and  wind  shear  (Troen  and  Mahrt,  1986).  The  formulation  is  based  on  bulk 
similarity  considerations.  The  depth  of  the  layer  is  represented  in  terms  of  modified  bulk 
Richardson  number.  Over  oceans,  the  surface  layer  and  the  turbulent  mixed  layer 
parameterizations  are  the  same  as  over  land  but,  unlike  over  land,  there  is  no  subsurface 
layer. 

The  COLA  model  planetary  boundary  layer  over  land  is  based  on  similar 
considerations.  It  also  consists  of  three  layers:  soil  layer,  surface  layer  and  turbulent  mix 
layer.  In  this  model  the  effects  of  the  soil  layer  and  the  vegetation  cover  is  treated  by 
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taking  into  accoun'  physiological  and  biological  processes  by  a  biosphere  model  (Sellers 
et  ai,  1986).  In  the  surface  layer  the  aerodynamic  resistances  for  momentum,  sensible  and 
latent  heat  transfer  between  the  earth’s  surface  and  the  atmosphere  are  based  on  the 
Monin-Obukhov  similarity  theory  (Miyakoda  and  Sirutis,  1987).  The  turbulent  mixed 
layer  which  varies  in  height  is  based  on  level  2  second-order  closure  model  of  Mellor  and 
Yamada  (1982).  Over  oceans,  the  surface  layer  and  the  turbulent  mixed  layer 
parameterizations  are  the  same  as  over  land.  There  is  no  subsurface  layer. 

The  precipitation  due  to  large-scale  supersaturation  and  subgrid  scale  moist 
convection  are  simulated  in  both  the  models.  The  precipitation  occurs  in  a  grid  box  if  the 
moist  air  is  cooled  such  that  relative  humidity  exceeds  100%.  However  the  precipitation 
occurring  in  a  higher  level  may  not  all  reach  the  earth’s  surface.  The  falling  precipitation 
may  evaporate  if  the  lower  model  layers  are  dry.  The  deep  moist  convective  precipitation 
occurs  if  there  is  a  moisture  flux  convergence  and  the  atmospheric  column  is  conditionally 
unstable.  The  shallow  convection  predominantly  occurs  in  undisturbed  flow  in  absence  of 
large  scale  moisture  flux  convergence.  The  trade  wind  cummuli  under  a  subsidence 
inversion,  daytime  convection  over  land  are  typical  examples  of  shallow  convection.  The 
shallow  convection  transfers  heat,  momentum  and  moisture  from  boundary  layer  to  the  free 
atmosphere.  For  deep  moist  convective  precipitation  the  AFGL  model  uses  Kuo’s  scheme 
(Kuo  1965,  1974)  modified  by  Krishnamurti  et  al.  (1976)  whereas  the  COLA  model  uses 
Kuo’s  scheme  modified  by  .\nthes  (1977).  For  shallow  convection  the  AFGL  model  uses  a 
semi-empirical  formula  derived  from  GATE  data  (Mahrt  et  ai,  1987).  The  COLA  model 
implements  the  scheme  of  Tiedtke  (1983)  which  enhances  the  vertical  diffusion  of  heat  and 
momentum  if  the  levels  near  the  surface  are  conditionally  unstable. 

Both  models  include  interaction  between  radiation  and  model-generated  clouds.  A 
broad  band  emissivity  approach  in  solar  and  longwave  radiation  in  presence  of  clouds  by 
Liao  and  Ou  (1981)  is  implemented  in  the  AFGL  model.  According  to  their  altitude,  high, 
middle,  and  low  clouds  can  form  if  the  model  generated  relative  humidity  exceeds  a 
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pre-selected  critical  value  as  suggested  by  Geleyn  (1981).  Middle  and  low  clouds  are 
considered  optically  black  whereas  high  clouds  are  half  black.  Convective  clouds  and  their 
interaction  with  radiation  are  not  included  in  the  AFGL  model.  The  COLiA  model 
implements  the  longwave  radiation  scheme  of  Harshvardhan  et  al.  (1987)  and  shortwave 
radiation  scheme  of  Lacis  and  Hanson  (1974).  Two  types  of  clouds,  convective  clouds  and 
supersaturation  clouds,  are  generated  in  the  model  following  a  procedure  similar  to  that  of 
Slingo  (1980,  1987)  (see,  Hou,  1990).  The  supersaturation  clouds  divided  into  three  types, 
high,  middle  and  low  according  to  their  altitude.  The  supersaturation  cloud  amounts  are 
determined  from  the  model  generated  relative  humidity  and  vertical  velocity.  Deep 
convective  cloud  amount  is  determined  from  the  convective  rainfall  rate  predicted  by  the 
Kuo  (1974)  convective  parameterization.  The  optical  properties  of  clouds  are  determined 
from  liquid  water  content  and  cloud  temperature. 

In  the  AFGL  model  the  following  quantities  are  prescribed;  the  sea  surface 
temperature  (SST)  from  FGGE  data  for  the  month  of  January,  roughness  length,  surface 
albedo,  estimated  snow  depth  and  soil  moisture  content  for  the  first  and  second  layer. 
Surface  specific  humidity  is  set  equal  to  zero  in  this  study.  Canopy  wetness  is  also  set 
equal  to  zero.  Soil  temperature  at  3  meter  deep  is  constant  in  time  but  varies  from  place 
to  place. 

For  the  COLA  model  the  SST  is  prescribed  from  NMC  data.  The  biosphere  model 
computes  roughness  and  surface  albedo.  But,  type  of  vegetation  cover  is  prescribed  from 
climatology  for  the  month.  Soil  moisture  and  snow  cover  are  prescribed  initially  but  they 
are  predicted  in  the  model.  Sea  ice  is  prescribed  from  NMC  analyses. 

There  is  also  a  difference  in  the  manner  in  which  SST  is  prescribed  in  the  two  models 
In  the  AFGL  model  SST  is  prescribed  at  all  ocean  grid  points  without  any  changes  In  the 
COLA  model  the  SST  is  modified  depending  on  the  height  of  the  a  =  1  surface  over  the 
oceans,  a  =  \  level  in  the  model  is  determined  by  the  prescribed  orography  in  the  spectral 
form.  However  due  to  the  Gibbs  oscillations  there  are  non  zero  height  values  for  a  =  1 
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level  over  the  oceans.  These  non  zero  values  could  be  as  much  as  1/2  km  above  or  below 
the  sea  level.  The  SST  values  over  ocean  are  modified  assuming  a  uniform  6.5°K/km  lapse 
rate  everywhere.  For  example,  if  the  height  of  a  =  1  level  over  ocean  is  ±100m  the 
prescribed  SST  is  modified  by  adding  ?0.65°K  to  SST  analysis  form  observations.  This 
modification  improves  the  simulation  of  heat  and  moisture  fluxes  at  the  ocean  surface. 

3.  EXPERIMENTS  AND  ANALYSIS 

VVe  have  made  nine  medium-range  forecasts  with  both  the  AFGL  and  COLA  models. 
NMC  analysis  was  used  for  initialization  and  forecast  verification.  The  nine  dates  were 
chosen  from  January  1990  data.  The  sample  size  of  nine  appears  reasonable  for  the  study 
of  systematic  errors  for  medium-range  forecast  (see  i.e.,  Harr  et  al,  1983).  Identical 
nonlinear  normal  mode  initialization  was  used  for  the  two  models.  As  mentioned  in  the 
Introduction  the  purpose  of  the  study  is  to  evaluate  the  performances  of  the  two  models  in 
comparison  with  observation  and  identify  the  causes  of  deficiencies  by  comparing  the  errors 
in  the  two  models.  Here  we  have  estimated  the  systematic  errors  in  the  two  models  and 
also  estimated  the  statistical  significance  of  the  systematic  errors  so  that  we  consider  the 
significant  error  for  identifying  the  causes. 

Let  Xf  and  xr,  represent  the  forecast  and  analysis  field  at  the  verifying  time 
respectively.  Hence  the  error  field  of  the  variable  is 

X  =  Xf  -  X  ,.  (1) 

If  we  have  m  forerasts  corresponding  to  m  initial  values  coming  from  fairly  uniform 
climatic  conditions  The  ensemble  error  field  is,  x  =  :^  Ex  where,  m  =  9. 

X  is  then  the  ensemble  mean  of  error  or  the  systematic  error.  The  forecast  field  x  for 
any  starting  date  can  be  expressed  as, 

X  =  X  +  x'  (2) 

where  x'  is  the  transient  error.  The  total  mean  square  error, 

ISx^  =  x2  + ls(x')2.  (3) 
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is  the  sum  of  square  of  the  systematic  error  and  the  mean  square  transient  error.  The  zonal 
average  of  the  error  field  is, 


X/,  =  5^/''xdA  (4) 

Jo 

where  A  is  the  longitude.  An  average  of  over  a  latitude  belt  between  01  and  02  is, 
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In  order  to  determine  the  significance  of  the  systematic  error  we  have  computed  the 
students  "t"  statistic  assuming  that  the  forecast  and  observed  variables  are  independent 
and  normally  distributed  with  the  same  population  variance  as: 


^v+Vo 


.Jm-1 


(5) 


where,  v  =  ^  E(xf-Xf)2,  Vo  =  ^  E(xo-Xo)2 


and  Xo  =  51xo,  mean  of  the  analysis. 

The  t-values  can  be  tested  with  a  pre-selected  level  of  significance.  The  t-value 
computed  has  (2m-2)  or  16  degrees  of  freedom  assuming  that  the  forecast  error  is  normally 
distributed  at  each  grid  point  and  the  nine  cases  selected  are  statistically  independent.  If 
the  t-values  lie  outside  the  range  ±2*1  would  imply  that  the  systematic  error  is 
significantly  different  from  zero  at  a  grid  point  for  a  two-tailed  test  at  the  95%  level. 
However,  one  has  to  be  careful  in  interpreting  the  significance  of  this  test.  First,  all  nine 
cases  may  not  be  statistically  independent,  second  this  test  is  valid  only  at  a  grid  point  and 
not  for  a  field.  Because  x  is  highly  correlated  with  neighboring  grid  points  more 
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sophisticated  tests  have  to  be  employed  for  field  significance  (Livezy  and  Chen,  1983).  The 
"t"  test  used  here  may  still  be  used  as  a  guidance  rather  than  an  accurate  measure  of 
statistical  field  signiGcance. 

4.  RESULTS 

We  shall  begin  with  a  comparison  of  root  mean  square  errors  in  the  500  mb 
geopotential  height  forecasts  with  the  AFGL  and  COLA  models.  Figure  1  shows 
systematic,  transient  and  total  errors  in  day-1  through  day-10  forecasts  averaged  over 
22°N-86°N  region  for  the  two  models.  The  error  levels  and  error  growth  rates  of 
systematic,  transient  and  total  errors  in  the  two  models  are  strikingly  similar.  The  error 
growth  rate  for  the  first  seven  days  is  larger  than  that  for  the  last  three-day  period.  Most 
likely  the  error  in  some  spatial  scales  is  saturated  by  day  7.  The  systematic  error  which  is 
less  than  20m  for  day  1  forecast  grows  to  about  70m  by  day  10.  The  systematic  error  is 
almost  half  as  much  as  the  total  error.  This  provides  significant  room  for  improving  the 
forecast  if  the  systematic  errors  are  entirely  due  to  the  approximations  in  simulating 
dynamical  and  physical  processes.  The  errors  in  500  mb  geopotential  height  for  day  1  to 
day  10  forecasts  averaged  over  the  tropical  belt,  22°N-22°S  are  shown  in  Figure  2.  Here 
there  is  a  significant  difference  between  the  systematic  errors  in  the  two  models.  The 
systematic  errors  for  days  1  to  10  in  the  AFGL  model  is  almost  twice  as  much  as  in  the 
COL.A  model.  For  day  10  the  systematic  error  is  about  33m  in  the  COLA  model  whereas 
in  the  AFGL  model  it  is  close  to  68m.  The  level  of  systematic  errors  for  days  1  to  10  in 
tropics  is  about  the  same  as  in  extra-tropics  in  the  AFGL  model.  The  transient  error  is 
about  the  same  in  the  two  models.  In  both  models,  the  systematic  error  dominates  the 
total  error,  in  contrast  to  extra-tropics.  Similar  characteristics  were  noted  in  ECMWF 
forecast  model  (Heckley,  1985).  In  the  tropics,  error  grows  rapidly  and  saturates  sooner 
than  in  extratropics  (Shukla,  1981).  The  error  growth  rate  for  the  first  three  days  is  larger 
than  that  for  the  remaining  part  of  the  forecast  period. 
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500  mb  height  error  (22N  -  86N  ave) 


day 

cola  sys  err 

afal  sys  err 

afal  trans  err 
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^  cola  total  err 

—  —  -  afgl  total  err 

Figure  1  The  time  evolution  of  500  mb  geopotential  height  errors  for  AFGL  and  COLA 
models  averaged  over  extratropics  (22°  N-86'  N)  for  ten  days  forecasts. 

Units:  m. 
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500  mb  height  error  (22N  -  22S  ave) 


Figure  2.  Same  as  Figure  1  except  the  errors  are  averaged  over  the  tropical  belt 
(22”  N-22"  S). 
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The  geographical  structure  of  systematic  error  of  the  500  mb  geopotential  height  for 
day  10  forecast  is  presented  in  Figure  3.  The  two  models  have  similar  error  structures  in 
that  they  are  negative  in  the  tropics  and  positive  in  the  extratropics  except  over  Arctic 
ocean  north  of  Eurasia.  These  general  characteristics  are  also  common  to  other  models: 
NMC  (Kanamitsu  et  al,  1990),  GFDL  (Sirutis  and  Miyakoda,  1990)  and  ECMWF 
(Heckley,  1985).  The  error  character  is  dominated  by  zonal  wave  number  four  in 
midlatitudes  and  zonal  wave  number  one  in  high  latitudes.  The  phase  and  amplitude  of 
the  error  fields  of  the  two  models  are  similar  with  minor  exceptions.  The  magnitude  of 
negative  errors  around  45°  N  and  45°  S  in  the  COLA  model  are  larger  than  those  for  the 
AFGL  model.  Large  tropical  errors  in  the  AFGL  model  seen  in  Figure  2  are  apparent  in 
Figure  3.  An  80m  isopleth  covers  a  large  monsoon  region  over  Indonesia,  Borneo, 
Philippines,  Southern  Brazil  and  the  Atlantic  Ocean.  A  60m  isopleth  covers  a  very  large 
part  of  the  tropics.  Whereas  in  the  COLA  model  40m  isopleth  covers  a  large  region  of 
tropics.  The  errors  in  tropics  can  be  a  serious  problem  for  the  predictability  of  extratropics 
flow  beyond  10  days  because  the  errors  from  tropics  propagate  and  contaminate  the 
predictions  in  extratropics  (see.  Palmer  et  ai,  1990).  The  geographical  distribution  of  the 
transient  error  in  day  10  for  the  two  models  is  shown  in  Figure  4.  The  transient  error 
patterns  in  the  two  models  are  very  similar.  The  error  level  in  extratropics  is  an  order  of 
magnitude  larger  than  in  tropics.  Comparing  the  extratropical  transient  error  pattern  with 
that  of  systematic  error  in  Figure  3  we  find  the  level  of  transient  error  is  large  where 
absolute  magnitude  of  systematic  error  is  large.  This  characteristic  is  also  common  to 
other  models,  ECMWF  and  GFDL  (see,  Hollingsworth  et  al.,  1980).  The  total  error  for 
day  10  is  shown  in  Figure  5.  It  reflects  the  combination  of  error  structures  seen  in  Figures 
3  and  4. 

Finally  the  distribution  of  t-values  for  day  10,  which  is  a  measure  of  significance  of 
the  systematic  error  is,  shown  in  Figure  6.  The  systematic  error  may  be  considered 
significantly  different  from  zero  in  a  region  where  absolute  value  of  t  exceeds  2.1.  As  noted 
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The  systematic  errors  in  500  mb  geopotential  height  for  day  10  forecast 
Upper  panel  is  for  AFGL  model  and  lower  panel  is  for  COLA  model. 
Contour  interval  is  20  m. 


earlier,  this  is  not  an  accurate  test  of  Geld  signiGcance,  therefore  we  should  interpret  the 
results  with  some  caution.  Both  models  show  similar  characteristics,  that  is,  the 
systematic  error  in  the  tropics  is  signiGcantly  different  from  zero  whereas  for  large  regions 
of  extratropics  they  are  not.  There  are  some  regions  over  North  Atlantic  ocean  and  Siberia 
where  the  systematic  errors  are  probably  different  from  zero.  In  tropics,  the  level  of 
signiGcance  is  higher  for  the  AFGL  model  compared  to  that  for  the  COLA  model. 

Next  we  shall  examine  global  error  statistics  for  850  mb  temperature  averaged  over 
22  N-HG'  N  region  shown  in  Figiin*  7.  1  he  general  characteristics  of  systematic,  transient 
and  total  errors  for  the  two  models  are  very  similar  to  those  for  the  500  mb  geopotential 
height  in  Figure  1  The  levels  and  growth  rate  of  errors  in  the  two  models  are  about  the 
same  The  level  of  systematic  error  in  the  two  models  is  about  1.5°  K  for  day  1  and  by  day 
10  it  is  about  3.4'  K  The  systematic  error  for  the  Grst  two  days  of  the  forecast  is  larger 
than  the  transient  error.  As  in  the  case  of  500  mb  geopotential  height  error,  the  rate  of 
error  growth  in  850  mb  temperature  for  the  Grst  seven  days  is  larger  than  that  for  the 
remaining  part  of  the  forecast  period. 

The  850  mb  temperature  error  statistics  for  the  tropical  belt  between  22' N  and  22°  S 
are  shown  in  Figure  8.  The  transient  error  level  in  the  two  models  are  about  the  same, 
slightly  lower  than  I  K  for  day  1  and  slightly  higher  than  I'K  for  10  day.  The  systematic 
errors  in  the  two  models  differ  signiGcantly.  For  day  1  the  systematic  error  it  is  about 
2.8’ K  in  the  COL.\  no  . del  and  3.4°  K  in  the  AFGL  model.  For  day  10  it  is  about  3.8°  K  in 
the  COLA  model  whereas  it  is  about  5.2' K  in  the  AFGL  model.  The  error  growth  rate  of 
systematic  error  is  large  for  the  Grst  tw'o  days  of  forecast.  By  day  4  it  appears  that  the 
error  is  probably  saturated  in  most  of  the  spatial  scales.  The  error  characteristics  in 
850  mb  temperature  are  similar  to  those  in  500  mb  geopotential  height  Geld.  This  is  more 
dearly  seen  in  Figure  9  whirh  shows  the  geographical  pattern  of  the  systematic  errors  for 
the  two  models  There  are  negative  errors  in  the  tropics  and  positive  errors  in  extratropics 
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850  mb  ten^.  err  (22N  -  86N  ave) 
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The  time  evolution  of  850  mb  temperature  error  for  AFGL  and  COLA 
models  averaged  over  extratropics,  (22°  N-86' N).  Units:  °K. 
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The  systematic  errors  in  850  mb  temperature  for  day  10  forecast.  Upper 
panel  is  for  AFGL  model  and  lower  panel  is  for  COLA  model.  Contour 
interval  is  2°K. 


except  over  Arctic  Ocean  north  of  Eurasia.  However  there  is  a  difference  in  the  error 
structure  between  850  mb  temperature  and  500  mb  geopotential  height,  that  is,  the 
magnitude  of  error  in  tropics  and  extratropics.  The  magnitude  of  850  mb  temperature 
error  in  tropics  is  about  the  same  as  in  extratropics  except  in  Northern  Territories  of 
Canada  which  differ  from  that  of  500  mb  height  field.  A  4°K  isotherm  covers  a  large 
region  of  tropical  oceans  in  both  the  models  but  in  the  AFGL  model  forecast  there  are  large 
regions  of  6°  K  isotherms  over  the  Pacific  and  Indian  Oceans  and  a  small  region  over  the 
Atlantic  Ocean. 

Similarities  in  Figure  3,  the  500  mb  geopotential  systematic  error  and  Figure  9,  the 
850  mb  temperature  systematic  error,  are  expected  because  of  the  hydrostatic  relation. 
The  general  pattern  of  these  two  figures,  that  is,  cold  in  tropics  and  warm  in  extratropics 
suggest  that  the  baroclinic  properties  of  the  models  are  going  to  differ  from  that  of  the 
atmosphere  because  of  weaker  north-south  temperature  gradient.  Some  of  the  error 
structure  in  the  extratropics  in  Figure  3  may  be  a  consequence  of  growth  rate  and 
propagation  speed  of  some  waves  in  the  model  differing  from  those  of  the  atmosphere 
(Lambert  and  Merilees  1978).  Again  this  characteristic  is  also  common  to  other  forecast 
models  and  GCMs  (Kanamitsu  et  al  1990;  Heckley  1985;  Sirutis  and  Miyakoda  1990).  The 
significance  of  850  mb  temperature  systematic  error  pattern  is  shown  in  Figure  10.  The 
errors  in  tropics  appear  highly  significant  in  both  the  models.  There  are  also  some  regions 
in  extratropics  such  as  over  Siberia,  Greenland  and  North  Atlantic  Ocean  where  the  errors 
are  probably  different  from  zero. 

The  errors  in  the  850  mb  relative  humidity  averaged  over  the  globe  (86°N-86°S)  for 
the  two  models  is  shown  in  Figure  11.  The  systematic  error  in  the  COLA  model  is  slightly 
lower  than  that  in  the  AFGL  model  for  the  two  days  forecast  but  then  it  increases  and  by 
day  10  it  is  about  2%  higher  than  that  in  the  AFGL  model.  The  transient  error  in  AFGL 
model  is  larger  than  that  in  the  COLA  model  and  as  a  result  the  total  error  in  the  AFGL 
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Figure  11. 


The  time  evolution  of  850  mb  relative  humidity  errors  in  per  cent  for 
AFGL  and  COLA  models  averaged  over  the  globe  (86°S-86'N). 
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Figure  12.  The  systematic  errors  in  850  mb  relative  humidity  for  day  10  forecast. 

Upper  panel  is  for  AFGL  model  and  the  lower  panel  is  for  COLA  model. 
Contour  interval  is  10%. 
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Figure  13. 


Same  as  Figure  12  except  for  t-values.  Contour  interval  is  2. 


model  is  larger  than  that  in  the  COLA  model.  The  error  growth  rate  is  much  higher  for 
the  first  four  days  of  forecast  compared  to  the  remaining  period  which  is  characteristic  of 
error  growth  in  the  tropics.  Geographical  distributions  of  the  systematic  errors  for  day  10 
in  the  two  models  are  shown  in  Figure  12.  In  the  AFGL  model,  negative  errors  are  largely 
over  the  oceans.  The  regions  of  negative  errors  are  over  the  north  and  southeastern  part  of 
the  Pacific  Ocean,  south  and  north  of  the  eastern  part  of  Atlantic  Ocean,  and  equatorial 
Indian  Ocean.  Over  land,  errors  are  largely  positive  except  for  small  regions  over  eastern 
Europe,  northeastern  Siberia  and  Northern  Territories  in  Canada.  In  the  COLA  model,  the 
errors  are  mostly  positive  except  small  regions  over  the  Pacific  Ocean  west  of  California, 
northwestern  Brazil,  eastern  Soviet  Union,  Saudi  Arabia  and  in  the  Indian  Ocean  south  of 
Indonesia.  If  the  relative  humidity  errors  are  largely  due  to  errors  in  the  850  mb 
temperature  one  would  expect  mostly  positive  errors  in  the  tropics  and  negative  errors  in 
the  extratropics  because  the  saturation  vapor  pressure  decreases  with  decreasing 
temperature.  In  the  COLA  model  the  temperature  errors  may  be  responsible  for  the 
relative  humidity  errors  to  some  extent.  The  errors  in  the  AFGL  model  may  be  caused  by 
the  fact  that  the  vertical  flux  of  moisture  over  oceans  is  not  properly  simulated  by  the 
boundary  layer  physics  or  by  the  shallow  convection.  A  comparison  of  tropical  relative 
humidity  errors  in  the  two  models  suggests  that  the  AFGL  model  is  drier  than  the  COLA 
model  whereas  the  former  is  colder  in  tropics  than  the  latter.  Hence  the  specific  humidity 
negative  errors  in  the  AFGL  model  are  probably  larger  because  the  saturation  vapor 
pressure  decreases  with  temperature. 

Errors  in  relative  humidity  or  specific  humidity  can  have  serious  consequences  on 
prediction  beyond  10  days  because  of  their  influence  on  other  physical  processes  such  as 
cloud-radiation  interaction  and  precipitation.  The  drier  atmosphere  (or  negative  specific 
humidity  error)  is  less  opaque  to  the  infrared  radiation  and  as  a  result  it  will  cool  further. 
Supersaturation  clouds  are  simulated  from  model  generated  relative  humidity.  Positive 
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relative  humidity  errors  are  likely  to  generate  more  clouds.  Assuming  that  the  albedo 
effect  (cooling)  of  clouds  dominates  the  greenhouse  effect  (warming),  more  clouds  will  cool 
the  atmosphere.  Whereas  the  negative  relative  humidity  errors  will  warm  the  atmosphere. 
The  moisture  convergence  and  conditional  instability  are  necessary  conditions  for 
convective  precipitation.  Precipitation  due  to  supersaturation  is  determined  from  the 
relative  humidity.  The  errors  in  relative  humidity  (specific  humidity)  can  affect  the 
amount  of  precipitation  and  hence  that  of  latent  heal  of  condensation. 

The  t-values  for  the  relative  humidity  are  shown  in  Figure  13.  The  area  of 
significant  systematic  errors  in  the  AFGL  model  is  less  than  that  of  the  COLA  model.  The 
magnitudes  of  systematic  errors  in  the  AFGL  model  are  slightly  lower  than  those  of  the 
COLA  model  but  at  the  same  time  the  transient  errors  in  the  AFGL  model  are  larger  than 
those  of  the  COLA  model  as  seen  in  Figure  11. 

The  errors  in  150  mb  temperature  averaged  over  extratropics  22°N-86°N,  are  shown 
in  Figure  14.  The  systematic  error  in  the  AFGL  model  is  about  the  same  as  that  of  the 
COLA  model  for  the  first  two  days  of  the  forecast  but  then  it  decreases  and  by  day  10  it  is 
about  0.4'’ K  less  than  that  of  the  COLA  model.  However  the  transient  error  in  the  AFGL 
model  is  larger  than  that  of  the  COLA  model.  The  error  growth  rate  characteristics  of 
150  mb  temperature  are  similar  to  that  of  850  mb  temperature.  The  growth  rate  is  larger 
for  the  first  seven  days  compared  to  that  of  the  remaining  period.  The  geographical 
patterns  of  the  systematic  error  in  the  two  models  are  shown  in  Figure  15.  Extratropical 
error  patterns  in  the  two  models  are  similar.  They  are  mostly  negative  except  over 
northern  Eurasia.  The  major  differences  are  in  the  tropical  error.  In  the  AFGL  model 
they  are  mostly  positive  and  large.  In  the  COLA  model  they  are  mostly  negative  with 
some  positive  error  regions  over  India,  Indian  Ocean,  south  Atlantic  Ocean  and 
southeastern  Pacific  Ocean.  The  significance  of  150  mb  temperature  error  is  shown  in 
Figure  16.  Tropical  systematic  errors  in  the  AFGL  model  appear  significantly  different 
from  zero  whereas  in  the  COLA  model  there  are  some  regions  where  the  error  is  not 
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Figure  14  The  time  evolution  of  150  mb  temperature  error  for  AFGL  and  COLA 

models  averaged  over  extratropics  (22°  N-86°  N).  Units:  °  K. 
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Figure  15. 


The  systematic  errors  in  150  mb  temperature  for  day  10  forecast.  Upper 
panel  is  for  AFGL  model  and  lower  panel  is  for  COLA  model.  Contour 
interval  is  2”  K. 
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Figure  16 


Same  as  Figure  15  except  for  t-values.  Contour  interval 


different  from  zero.  In  both  the  models  errors  in  the  southern  hemisphere  extratropics  are 
significant.  In  the  northern  hemisphere  extratropics  there  are  only  a  few  regions  where  the 
error  is  significant. 

Now  we  shall  consider  the  error  structure  and  their  statistical  significance  in  zonally 
averaged  temperature  and  zonal  winds.  The  systematic  errors  of  zonally  averaged 
temperature  for  the  AFGL  and  COLA  day  10  forecasts  are  shown  respectively  in 
Figures  17a  and  Fig.  17b.  Both  models  have  similar  error  structure  in  that  the  troposphere 
is  cold  except  in  high  latitudes  of  lower  troposphere  and  equatorial  lower  stratosphere.  The 
only  major  difference  in  the  two  figures  is  the  magnitude  of  tropical  stratospheric  error 
which  is  larger  in  the  AFGL  model  than  in  the  COLA  model.  The  characteristics  of  error 
structure  in  850  mb  temperature  in  Figure  9  that  cold  tropics  and  warm  extratropics  can 
be  seen  in  the  zonally  averaged  error  field.  As  mentioned  earlier,  such  error  structure  will 
reduce  the  baroclinicity  of  the  model  atmosphere.  The  vertical  error  structure,  cold  in 
lower  troposphere  and  warm  in  lower  stratosphere  will  reduce  the  convective  instability  of 
the  tropical  model  atmosphere.  However,  these  characteristics  are  not  only  common  to 
these  two  models  but  also  in  other  forecast  and  GCM  models  (Kanamitsu  et  a/.,  1990; 
Sirutis  and  Miyakoda,  1990;  and  Ileckley,  1985).  As  shown  by  the  distributions  of 
student’s  "t"  statistic  in  Figures  18a  and  18b,  these  systematic  errors  are  statistically 
significant. 

The  systematic  error  in  zonally  averaged  zonal  wind  for  the  AFGL  and  COLA 
models  is  shown  in  Figures  19a  and  19b  respectively.  The  error  structures  in  the  two 
models  are  similar  in  that  the  errors  have  a  barotropic  structure.  Both  models  have 
easterly  bias  in  tropics  and  westerly  bias  in  southern  hemispheric  subtropics  and  again 
easterly  bias  in  high  latitude  of  the  southern  hemisphere.  Both  models  have  positive  errors 
in  high  latitudes  of  the  northern  hemisphere  and  negative  errors  around  45°  N.  The  COLA 
model  has  westerly  bias  in  subtropics  of  the  northern  hemisphere  but  in  the  AFGL  model 
they  are  only  in  the  stratosphere  and  in  the  troposphere  it  has  an  easterly  bias.  Prior  to 
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Figure  18b. 


Same  as  Figure  17b  except  for  t-values.  Contour  interval  is  2. 
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the  inclusion  of  explicit  parameterization  of  gravity  wave  drag  the  magnitude  of  these 
errors  were  much  larger  (see,  Kirtman,  et  al  1991).  The  error  structure  of  the  zonal  winds 
are  consistent  with  errors  in  the  zonally  averaged  temperature  in  that  they  approximately 
satisfy  the  thermal  wind  equation.  Finally,  the  structure  of  the  t-values  for  the  AFGL  and 
COLA  models  are  respectively  shown  in  Figures  20a  and  20b.  The  major  feature  of  the 
error  structure  is  significant.  However  there  are  regions  in  the  vicinity  of  zero  isopleth  the 
errors  are  not  statistically  significant. 

5.  SUMMARY  AND  CONCLUSIONS 

We  have  made  nine  ten-day  forecasts  with  the  two  models,  AFGL  and  COLA,  to 
study  the  systematic  errors  with  the  intent  of  identifying  the  sources  of  these  errors.  Nine 
initial  dates  were  chosen  from  January  1990  such  that  they  are  synoptically  independent. 
The  geographical  patterns  of  systematic  errors  and  their  statistical  significance  are 
computed  for  500  mb  geopotential  height  field,  850  mb  temperature  and  relative  humidity 
and  150  mb  temperature.  Additionally,  vertical  and  meridional  distributions  of  systematic 
errors  and  their  statistical  significance  were  computed  for  zonally  averaged  temperature 
and  zonal  wind.  The  systematic  errors  in  the  two  models  have  sonse  common  features.  For 
example,  the  errors  in  500  mb  geopotential  height  are  negative  in  tropics  and  positive  in 
extratropics.  Consistent  with  500  mb  height  errors,  the  850  mb  temperatures  are  cold  in 
tropics  and  warm  in  extratropics  compared  to  NMC  analysis  The  zonally  averaged 
temperatures  in  tropics  are  cold  in  the  troposphere  and  warm  in  the  lower  stratosphere.  In 
the  extratropics  the  temperatures  are  warm  in  lower  troposphere  and  cold  in  the  upper 
troposphere  compared  to  the  observations.  The  systematic  errors  of  zonally  averaged  zonal 
winds  have  a  barotropic  character.  By  and  large,  the  errors  are  negative  (easterly  bias)  in 
the  tropics,  middle  latitudes  and  high  latitudes  of  the  southern  hemisphere.  The  positive 
errors  (westerly  bias)  appear  in  the  subtropics  and  high  latitudes  of  the  northern 
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Figure  19a. 


The  systematic  errors  in  zonally  averaged  zonal  wind  for  day  10  forecast 
with  AFGL  model.  Contour  interval  is  2  m/sec. 
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hemisphere.  The  errors  in  the  zonally  averaged  zonal  wind  are  consistent  with  those  of 
temperature  errors  to  approximately  satisfy  the  thermal  wind  relation.  Extratropical 
systematic  errors  in  day  10  forecast  of  500  mb  geopotential  height  and  850  mb 
temperatures  in  both  the  models  are  larger  than  those  in  tropics  but  they  are  not  all 
significantly  different  from  zero  because  the  transient  errors  are  also  large.  The  tropical 
systematic  errors  in  both  the  models  appear  significantly  different  from  zero.  One  of  the 
differences  between  the  error  characteristics  of  the  two  models  is  the  magnitude  of  error  in 
the  tropics.  The  AFGL  model  tropical  error  magnitudes  in  500  mb  geopotential  height, 
and  850  mb  and  150  mb  temperature  are  larger  than  the  corresponding  errors  in  the  COLA 
model  in  some  cases  twice  as  large.  There  are  large  differences  in  the  systematic  error 
structure  in  850  mb  relative  humidity  between  the  two  models.  In  the  AFGL  model 
negative  errors  are  largely  over  ocean  and  positive  errors  are  over  land.  In  the  COLA 
model  the  errors  are  largely  positive  with  some  small  regions  of  negative  errors.  In  the 
AFGL  model  the  relative  humidity  errors  appear  rather  serious.  The  model’s  tropical 
tropospheric  temperature  has  a  cold  bias.  Because  the  saturation  vapor  pressure  decreases 
with  temperature,  one  would  expect  positive  error  in  the  relative  humidity  as  in  the  COLA 
model.  The  negative  errors  imply  that  lower  troposphere  over  tropical  oceans  is  very  dry. 

There  are  two  major  differences  in  the  manner  in  which  the  physics  is  treated  in 
these  models.  First  is  that  the  AFGL  model  does  not  include  radiation  interaction  with 
deep  convective  clouds.  Second  is  the  manner  in  which  the  SST  is  prescribed.  The  a  =  1 
surface  in  the  model  is  determined  by  the  prescribed  orography  in  the  spectral  form.  Due 
to  the  Gibbs  oscillations  there  are  some  non-zero  height  values  over  the  ocean.  The  SST 
prescribed  form  observations  create  fictitious  horizontal  temperature  gradients  on  a  =  1 
surface  because  it  is  lower  or  higher  than  the  sea  level.  This  can  have  an  adverse  effect  on 
heat  and  moisture  fluxes  between  the  surface  and  the  atmosphere.  To  minimize  this  effect 
the  SST  in  the  COLA  model  is  interpolated  to  a  =  1  surface  assuming  uniform  6.5°K/km 
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lapse  rate.  This  modification  in  SST  is  not  included  in  the  AFGL  model,  but  its  effect  is  to 
some  extent  minimized  by  the  smoothed  orography.  The  tropical  temperature  errors  in  the 
AFGL  could  be  reduced  by  including  the  convective  cloud-radiation  interaction.  Heating 
at  the  cloud  base  in  the  lower  troposphere  and  cooling  at  the  cloud  top  in  the  lower 
stratosphere  will  improve  the  temperature  forecast  in  the  tropics.  Interpolating  the 
prescribed  SST  to  a  =  1  level  may  improve  the  error  in  relative  humidity  and  temperature 
over  oceans.  The  AFGL  model  includes  a  shallow  convection  scheme  based  on  an  empirical 
relation  derived  from  the  GATE  data  (Mahrt  et  al  1987).  It  is  possible  that  this 
parameterization  is  not  as  effective  as  the  atmosphere  in  transporting  heat  and  moisture 
from  the  surface  to  the  atmosphere.  The  lack  of  vertical  transport  of  heat  and  moisture 
may  be  responsible  for  the  large  errors  in  850  mb  temperature  and  relative  humidity. 
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APPENDIX  A 


SILHOUETTE  OROGRAPHY  REPRESENTATION 


Highest  resolution  global  orography  data  available  at  the  present  is  the  10- 
minute  elevation  data  from  U.S.  Navy.  The  fine  orographic  structure  in  the  Navy 
data  cannot  be  explicitly  represented  in  models  of  present  horizontal  resolution.  As 
mentioned  in  the  Introduction  the  mean  orography  computed  as  an  arithmetic  aver¬ 
age  of  points  over  the  model  grid  underestimates  the  blocking  effects  of  high 
mountains.  In  order  to  achieve  the  proper  large  scale  effects  of  mountains  the  so 
called  silhouette  orography  is  computed  by  taking  an  average  of  maximum  peaks  of 
mountain  profiles  in  the  east-west  and  north-south  directions  over  a  grid  box,  that 
is. 
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where  N  and  M  are  the  number  of  subgrid  profiles  of  mountain  within  a  grid  box 
in  X  and  y  directions,  respectively.  ZI  is  maximum  height  of  the  profile  in  the  y-z 
plane,  and  similarly,  ZJ  is  the  maximum  height  of  the  profile  in  the  x-z  plane. 


The  silhouette  orography  computed  from  equation  (A-1)  is  then  represented 
in  terms  of  the  spherical  harmonics.  The  number  of  harmonics  used  in  the 
representation  depends  on  the  horizontal  resolution  of  the  model.  Rhomboidal 
truncation  at  wave  number  30  is  used  in  the  AFGL  model.  The  truncation  pro¬ 
duces  Gibbs  oscillations  in  the  representation.  The  amplitudes  of  the  oscillations 
over  ocean  could  be  as  large  as  I  km.  To  reduce  the  amplitude  of  these  oscillations 
a  scale  selective  smoother,  suggested  by  Hoskins  (1980),  is  applied  to  the  spectral 
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coefficients.  The  smoother  is, 
exp(-  k(  n(  n+  1 )  )2) 

where  k  is  selected  such  that  the  highest  retained  coefficient  is  reduced  to  0.1  of  its 
initial  value,  and  n  is  the  degree  of  associated  Legendre  function.  The  smoothing 
function  reduces  the  amplitude  of  oscillations  over  ocean  and  yields  smoother 
topography  field  over  continents  without  significantly  reducing  the  enhanced  effect 
of  silhouette  orography. 
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APPENDIX  B 


PARAMETERIZATION  OF  OROGRAPHIC  GRAVITY  WAVE  DRAG 

As  mentioned  earlier,  depending  on  atmospheric  stability  and  vertical  wind 
shear,  gravity  waves  can  propagate  vertically  to  great  heights  until  they  are 
absorbed  and/or  reflected  by  critical  layers,  or  they  become  unstable  as  a  conse¬ 
quence  of  convective  or  shear  instabilities,  in  which  case  the  gravity  waves  dissi¬ 
pate.  These  waves  transport  momentum  vertically  and  deposit  it  where  they  are 
dissipated.  The  level  where  they  are  dissipated  the  momentum  is  lost  by  the  large- 
scale  flow  and  it  is  brought  down  to  the  earth’s  surface  where  it  is  deposited.  The 
parameterization  of  these  effects  in  GCMs  and  NWP  models  are  based  on 
simplified  theoretical  concepts  and  observational  evidence.  The  parameterization 
consists  of  determining  the  drag  due  to  gravity  waves  at  the  surface  and  its  vertical 
variation  in  the  atmosphere.  There  are  two  approaches  to  parameterize  the  wave 
drag  at  the  surface.  One  is  based  on  linear  theory  (Palmer  et  a/., 1986;  McFar- 
lane.1987)  and  the  other  is  based  on  nonlinear  theory  (Pierrehumbert,  1987).  The 
vertical  variation  of  the  wave  drag  in  the  atmosphere  depends  on  critical  layers  and 
convective  or  shear  instabilities  (Palmer  et  ai.  1986;  Helfand  et  a/., 1987).  Here  we 
present  the  parameterizations  of  the  surface  drag  based  both  on  the  linear  and 
nonlinear  theories. 

1  GRAVITY  WAVE  DRAG  AT  THE  EARTH’S  SURFACE 
a)  Linear  Theory 

Let  us  consider  a  two  dimensional,  adiabatic,  inviscid  flow  over  a  small  moun¬ 
tain,  50  to  100  km  wide.  The  governing  equations  for  the  mountain  wave  can  be 
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written  as: 
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We  shall  express  each  variable  as  a  sum  of  a  basic  state  and  a  perturbation, 
where  the  basic  state  is  given  by, 
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Linearizing  equations  B-1  to  B-5  and  assuming  a  steady  state  condition,  we 


have. 
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(B-8) 


rr^p'  .  'dp  .  9w\ 


0-^+w'^  =  c2(g|el+w'^) 
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(B-9) 


Equation  (B-9)  may  be  rewritten  as, 

p  dx  p  dz  c^p  dx 


(B-10) 


where  the  last  term  can  be  neglected,  which  implies  that  the  pressure  variations  are 
not  important  in  generation  of  density  anomalies  for  low  frequency  motion.  Define 


P  =  _  _I_  =  J_d0 
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which  is  a  measure  of  the  static  stability.  The  equations  B-6  to  B-10  can  be  reduced 
to  obtain  a  single  equation  for  vertical  velocity  w'(x,z). 
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^  U  dz  w  dz 


(B-12) 


■vhere 


S= 

dz 


(B-13) 


Let 


2..,' 


w  =  [p  /  p,)  w 


and  neglect  the  last  term  of  B-12,  since  - 

becomes 


10^ 

300^ 


(B-14) 


<<  1.  Equation  B-12  then 
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Wxx+  Wz2+  l^(z)w  =  0 
where  Scorer  parameter  is  defined  as. 


(B-15) 


u2  u  4  2  ^  u 

Equation  B-15  is  well  known  as  the  linearized  steady  state  equation  for  the 
vertical  motion  of  gravity  waves.  In  practice,  l^(z)  is  usually  dominated  by  the  first 
term,  i.e.,  the  buoyancy  force  term.  Only  in  regions  of  strong  shear  will  the  term 
Ujj/U  become  important.  Neglecting  the  S  terms  is  equivalent  to  making  the 
Boussinesq  approximation,  that  is,  density  variations  are  only  important  as  they 
affect  the  buoyancy.  To  determine  the  lower  boundary  condition,  the  flow  is 
assumed  to  follow  the  terrain  at  the  ground,  which  is  described  by  sinusoidal  topog¬ 
raphy,  h*(x)=  h^sinkx.  The  streamline  slope  equals  the  terrain  slope: 

u>  u/'  dh* 

(atz=h*(x))  (B-16) 

u  U  +  u  dx 

=  h^kcoskx.  (B-17) 

Assuming  the  amplitude  of  topography  is  small  (as  is  disturbance  u'),  the 
lower  boundary  condition  may  be  expressed  as: 

w' =  w  =  Uh^kcoskx  (B-18) 

In  order  to  solve  the  equation  B-15,  we  shall  separate  the  variables,  that  is,  let 

w(x,z)  =  4>](z)  coskx-^  <I)2(z)  sin  kx.  (B-19) 

Substituting  B-19  into  B-15  we  get, 

0^3 -t-  ( l\z)  ~k^)0=  0.  (B-20) 

Both  Oi  and  <1>2  satisfy  equation  (B-20).  To  satisfy  the  lower  boundary  condition,  k 
should  be  the  same  as  the  terrain  wave  number.  The  sign  of  the  term  in 
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parentheses  of  the  equation  is  important  because  it  determines  whether  the  station¬ 
ary  gravity  waves  are  internal  or  external.  Here  we  consider  only  vertically  pro¬ 
pagating  internal  gravity  waves,  i.e.,  k^<  If  1  is  treated  as  piecewise  constant  as 
Scorer  did,  the  solution  of  (B-20)  is  obtained  as  follows: 

ji_  X 

<I>(z)  =  A  sin  (1^-  k^)  ^z+  Bcos(l^-  k^)  ^z 

combining  this  with  (B-19)  gives 

w(x,z)  =  C  cos  [  kx  (1^  -  k^)  ^  z  ]  +  D  cos  [  kx  -  (1^-  k^)  ^  z  ] 

1  1 

-t-  Esin[kx-t-  (1^-  k^)  2  z  ]  +  F  sin  [  kx  -  (l^-  k2)2z]  (B-21) 

Applying  the  lower  boundary  condition  implies 

E+  F  =  0 
C+  D  =  h^kU 

The  so-called  radiation  condition  is  imposed,  i.e.,  that  there  be  no  components  of 
the  flow  which  radiate  energy  downward. 

D  =  F  =  0 

This  last  condition  implies  that 
E=0 
and 

C  =  h^kiJ. 

Thus,  equation  (B-21)  becomes 

w(x,z)  =  U  hf„k  cos  f  kx  ■+ (l2- k^)  2  z  J.  (B-22) 

If  the  hydrostatic  assumption  is  employed,  it  can  be  shown  that  w^^  is  eliminated  in 
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(B-15),  such  that  disappears  in  (B-20)  leading  to  the  solution: 


w(x,z)  =  U  h,„kcos(kx+ Iz) 


(B-23) 


where. 


P.  2  P* 

w(x,z)  =  —  w(x,z)  =  —  U  hmkcos(kx+lz) 

IpJ  IPJ 


Under  the  Boussinesq  approximation,  the  continuity  equation  becomes 

1^1+ 1^1=  0. 

dx  dz 


(B-24) 


So  there  exists  stream  function  which  can  be  defined  by 

w'=U^ 

dx 


u'= 

8z  ■ 

Thus 


r 

y  =  J  dx  =  h  sin  (kx  +  Iz) 


(B-25) 


p.  2 


where,  h  =  J  h^,  is  the  displacement  amplitude  of  stream  function, 
dz 

=  -  U  hj,  sin  (kx  +  Iz)  -  Uj  h  sin  (kx  +  Iz) 


-  U  h  lcos(kx+  Iz) 


(B-26) 


The  upward  flux  of  momentum  averaged  over  one  horizontal  wave  length  L 


now  can  be  written  as 
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1^ 

fycTw'  =  -^1  u'w'< 


dx 


=  -  (h  Uk  (  Ujh  +  U  hj,) )  J  sin  (kx+  Iz)  cos(kx+  Iz)  dx 


_  £fT2h2i 


-U^h^klj  cos(kx+ lz)^dx 
^  0 


Since  J  sin(kx+  Iz)  cos(kx+  Iz)  dx  =  0, 


U 

pFiTw' =  - -p-O^h^kl  J  -^  ( 1  -  cos(2(kx+ Iz)))  dx 


=  -  -^-ij^h^kl 
2 


Recall, 


(B-27) 


.2  =  ll.  =  Ni 

where  N  is  the  Brunt  -  Vaisala  frequency.  Substituting  this  into  (B-27)  we  obtain 
the  gravity  wave  drag  at  the  surface  as  follow: 


T  =  -  p  u  w 

=  -S^kh^UN 

2 

=  Kph^iJ  N 


(B-28) 


where  K  =  —k.  The  horizontal  characteristic  subgrid-scale  wavelength  is 

L  =  ~  =  125  km 
k 


so  that 


K  =  2.5x  10  ^  m 


-5  „-l 
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h^,  the  variance  of  subgrid-scale  orography,  is  computed  from  10-minute  elevation 
data  from  the  U.  S.  Navy.  U,  p  and  N  are  the  model  variables.  In  order  to  avoid 
numerical  instabilities,  h  is  limited  to  a  maximum  value  of  400  m.  In  addition,  the 
linear  theory  requires  that  Froude  number, 

Fr  =  0.8 

U 

If  Fr  exceeds  0.8  nonlinear  effects  become  important.  Therefore,  wherever  Fr 
exceeds  0.8  the  surface  drag,  which  depends  on  h^,  is  reduced  by  a  factor  (0.8/ Fr)^. 
That  is. 


Max  (h)  =  min  {0.8  4(X)m}. 

N 


b)  Nonlinear  Theory 

The  parameterization  of  surface  wave  drag. 


X  =  xpNUh^ 


(B-29) 


based  on  the  linear  theory,  is  valid  only  for 


Fr  =  0.8. 


Pierrehumbert(  1987)  has  argued  that  very  close  to  the  surface  where  U  is  small  and 
in  the  vicinity  of  mountains  with  large  h  Fr  exceeds  0.8  and  hence  the  linear 
theory  is  not  valid.  He  proposed  a  parameterization  of  surface  wave  drag  to  include 
situations  where  Fr  is  greater  than  0.8.  It  is  based  on  dimensional  considerations 
and  results  from  numerical  experiments. 


Ixl  =  ou'w'  = 


j  pu'w'dx 


(B-30) 


where.  Ax  is  grid  length  of  the  model.  Equation  (B-30)  can  be  written  as 
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(B-31) 


lx  I  «  -^p(u1[w']L 
Ax 


where,  L  is  the  length  scale  of  mountain.  Let, 


[u'l  =  U 


[w']  = 


(B-32) 

(B-33) 

(B-34) 


Substituting  (B-32)  and  (B-34)  in  (B-31),  we  get 


Itl  « 


Ax  N 


or 


(B-35) 


'^'=^^G(Fr)  (B-36) 

Ax  N 

where,  G(Fr)  is  a  monotonically  increasing  function  of  Fr.  The  functional  form  of 
G(Fr),  determined  from  numerical  experiments,  is 


G(Fr)  = 


Fr^ 

1-h  Fr^ 

V,  J 


If  there  are  m  mountains  in  a  grid  box. 


Ixl 


m 

5^1 

Fr^ 

Ax 

l-l-  Fr^ 

k  J 

(B-37) 


• 

Fr^ 

L 

h  * 

l-t-  Fr^ 

V  J 

(B-38) 


where,  L  corresponds  to  wavelength  of  monochromatic  wave  in  the  direction  of  the 
surface  wind.  The  equation  (B-38)  is  valid  for  a  wide  range  of  values  of  Fr.  For 
values  of  Fr  less  than  0.8  where  linear  theory  is  valid,  the  equation  (B-38)  reduces 
to 
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It  I  =  pKUNh^ 

which  is  the  same  as  (B-28)  the  parameterization  based  on  linear  theory. 

Pierrehumbert(1987)  has  shown  that  the  nonlinear  effects  can  lead  to  large 
changes  in  the  character  of  the  flow  and  amplification  of  wave  drag  close  to  the  sur¬ 
face.  To  account  for  the  nonlinear  amplification  a  base  layer  is  defined  to  compute 
the  surface  wave  drag.  The  base  layer  is  approximately  the  first  third  of  the  atmo¬ 
sphere.  We  have  taken  this  layer  to  be  from  surface  to  642  mb.  The  surface  drag  is 
computed  from  equation  B-38  by  using  mass  weighted  model  variables  over  the  the 
depth  of  the  base  layer.  L  is  equal  to  125  km  as  before.  Also,  in  this  case  the  max¬ 
imum  value  of  h  is  restricted  to  400  m  for  numerical  stability. 

2.  GRAVITY  WAVE  DRAG  IN  THE  ATMOSPHERE 

Eliassen  and  Palm(1961)  have  shown  that  in  absence  of  transience  and  tur¬ 
bulent  dissipation  the  momentum  flux  is  independent  of  height.  In  such  a  case  the 
momentum  flux  is  deposited  in  the  top  layer  of  the  model  and  hence  there  is  no 
body  force  in  the  model  atmosphere  because  the  vertical  variation  of  momentum 

flux,  vanishes.  The  momentum  flux  can  vary  in  vertical  if  one  or  more  of  the 
do 

following  conditions  occur:  (i)  wave-modified  Richardson  number,  Ri^,  is  less 
then  the  critical  Richardson  number,  Ri,.,  (ii)  shear  instabilities  (Ri  <  1/4)  create 
turbulence,  (iii)  convective  instabilities  (Ri  <  0)  create  turbulence  and  (iv)  critical 
level  where  direction  of  the  wind  is  perpendicular  to  surface  wind  or  greater  than 
90  degrees. 

The  wave-modified  Richardson  number  can  be  expressed  in  terms  of  the 
Richardson  number  for  the  undisturbed  state  of  the  atmosphere.  The 
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Brunt  -  V'aisala  frequency  is 


2  =  _£_  £?I 


Bj-  dz 


(B-39) 


where  the  subscript  T  refers  to  the  sum  of  the  background  flow  and  wave  contribu¬ 


tion,  i.  e. 


Bf  —0+0 


—  +  —  ). 
0+  0'  dz  dz 


(B-40) 


Under  the  adiabatic  approximation  the  thermodynamic  equation  can  be  written  as 


dx  dz 


(B-41) 


So  0' can  be  written  in  terms  of  stream  function. 


(B-42) 


0-t^  0'  dz  dz  ^  dz^ 


(B-43) 


Neglecting  the  second  order  derivative  term  and  approximating 


0+  0'  0 


we  get. 


0  dz  dz 


(B-44) 


Referring  back  to  equation  B-25  that 


y  =  h  sin  (  kx  +  Iz) 


(B-45) 
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thus 


Nj^  =  ( 1  -  h  1  cos  ( kx  +  Iz )  +  sin  (  kx  +  Iz ) ) 

=  N^(  1-  ■^^cos(kx+  Iz)  -  4"“^ sin  (kx+lz)) 

U  2  g 

>  N2(1--?^)  =  Nj  (B-46) 

where  N  =  4-4^  which  is  calculated  from  the  mean  flow.  The  term  - 

0  dz  U 

can  be  regarded  as  the  influence  at  the  phase  of  the  (monochromatic)  wave  for 

which  the  decrease  in  the  local  static  stability  is  maximized.  Finally,  the  wave 

modified  Richardson  number  is 


^im  = 


-^=R.(1-^) 
9U  /  ‘  U 

^  dz  ^ 


(B-47) 


This  implies  that  for  small  values  of  U  wave  breaking  will  occur  even  for  moun¬ 
tains  of  relatively  moderate  heights.  The  critical  Froude  number,  Fr<.,  for  breaking 


can  be  written  as 


Frc=  Nh^/U  =  1-  Ric/Ri=  1-  l/4Ri  (B-48) 

Once  the  gravity  wave  begins  to  break,  the  wave  saturation  hypothesis  of 
Lindzen  (1981)  is  invoked,  in  which  it  is  assumed  that  the  wave-induced  turbulent 
dissipation  prevents  the  displacement  amplitude,  h,  from  exceeding  its  critical  value 
hp,  that  is, 

\=  Fr^U/N  =  (Fr^U/Nh)  =  (Fr^/ Fr)h 
Replacing  h  by  h^  in  equation  (B-28)  we  have 

Xc  =  pUNKh/  =  pUNKh2(Frc/  Fr)^  =  t  (Fr^/  Fr)^  (B-49) 

This  implies  that  if  the  computed  value  of  Fr  at  a  given  level  exceeds  critical  value, 
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Fr,;, 

it  is  assumed  that  the  turbulence  generated  by  the  breaking  wave  will  reduce  the 
wave  momentum  flux  by  the  factor  (Frc/Fr)^.  However,  h  the  displacement  of 
isentropic  surface  on  a  particular  level  is  not  known.  The  following  procedure  is 
used  to  calculate  Fr.  From  equation  (B-28),  we  have 

X  =  KpCfNh^ 

h2=-4_=-^- 

KpUN  KpLTFr 

h  =  - - 

Kp  U^Fr 

Substitute  this  relation  of  h  into  the  definition  of  Froude  number  Fr=hN/U  we 
obtain  the  expression  for  the  square  of  the  Froude  number, 

Fr^=  — (B-50) 
KpU^ 

Very  close  to  the  surface  convective  and  shear  instabilities  occur  frequently 
and  dissipate  the  momentum.  In  order  to  allow  gravity  waves  to  propagate  verti¬ 
cally  through  boundary  layer  the  parameterization  based  on  nonlinear  theory  where 
the  base  layer  extends  from  surface  to  642  mb  the  momentum  flux  in  the  base 
layer  can  change  only  due  to  wave  saturation,  that  is,  if  Rin,<  Ric-  Here  the 
effects  of  shear  and  convective  instabilities  or  critical  levels  are  ignored. 

After  the  vertical  profile  of  the  momentum  flux  has  been  determined,  the 
influence  of  the  gravity  wave  drag  on  the  large  scale  flow  is  determined  by  the  vert¬ 
ical  wave  momentum  flux  divergence. 
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APPENDIX  C 


COMPUTER  CODE  FOR  OROGRAPHIC  GRAVITY 
WAVE  DRAG  PARAMETERIZATION 


SUBROUTINE  GWDD(PSFC,U,V,TT,CHUG,CHVG,IMX,KMAX,LATCO) 
DIMENSION  XTENS(96,19),YTENS(9«,I9) 

DIMENSION  RO(96,18),PP(96,l«),TENSlO(96,19) 

DIMENSION  SLEV(19),SLAY(18) 

DIMENSION  DZ(96,I8),PPP(96,19) 

DIMENSION  DRAGSF(96),XDRAG(96),YDRAG(96) 

DIMENSION  TT(IMX,KMAX),T(98,18). 

1  U(IMX,KMAX),V(IMX,KMAX),PSFC(IMX) 

DIMENSION  CHUG(IMX,KMAX),CHVG(IMX,KMAX) 

COMMON  /SF/  VAR(96,80),DOIT,UTEN(96,80,18), 

1  VTEN(96,80,18),TAUSX(96,80),TAUSY(96,80) 

INTEGER  DOIT 

REAL  LSTAR,NBAR,BV(18) 


C 

C  INPUT  PARAMETERS 

C 

C  PSFC  SURFACE  PRESSURE 

C  U  ZONAL  WIND 

C  V  MERIDIONAL  WIND 

C  TT  VIRTUAL  TEMPERATURE 

C  IMX  NUMBER  OF  GRID  POINTS  ALONG  LONGITUDE 
C  KMAX  LEVELS  IN  THE  VERTICAL 
C  LATCO  NUMBER  OF  GAUSSIAN  COLATITUDES 
C  VAR  OROGRAPHICAL  VARIANCE 

C 
C 

C  OUTPUT  PARAMETERS 

C 

C  UTEN  U  TENDENCY 

C  VTEN  V  TENDENCY 

C 
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MLONG=  96 

DO  123  1=  MLONG+  l.IMX 
DO  123  K=  l.KMAX 
CHUG(1,K)=  0.0 
CHVG(1.K)=  0.0 

123  CONTINGK 

DO  124  1=  l.MLONG 
DO  124  K=  l.KMAX 
T(I,K)=  TTd.K) 

124  CON  TIN  UK 


DO  2000  1=  l.MLONG 

PSFC(I)=  KXP(PSFC(I))*I0. 
2000  CONTINUK 


MLAT=  80 
NIAY=  18 
NLAYM1=  NLAY-1 
NLAYPU  NLAY+  1 
DLON=  360.0/96.0 
DLAT=  180.0/80.0 
GRAV=  9.81 

agrav=  1.0/ gray 

RGAS=  287.0 
PTOP=  10.0 
ICOUNT=  0 

SLEV(19)=  0.0 
SLEV(18)=  0.05 
SLEV(17)=  0.1 
SLEV(16)=  0.15 
SLEV(15)=  0.2 
SLEV(14)=  0.25 
SLEV(13)=  0.3 
SLEV(12)=  0.35 
SLEV(11)=  0.4 
SLEV(10)=  0.45 
SLEV(9)=  0.546 
SLEV(8)=  0.642 
SLEV(7)=  0.735 
SLEV(6)=  0.82 


SLEV(5)=  0.893 
SLEV(4)=  0.948 
SLEV(3)=  0.973 
SLEV(2)=  0.99 
SLEV(1)=  1.0 

SLAY(18)=  0.021 
SLAY(17)=  0.074 
SLAY(16)=  0.124 
SLAY(15)=  0.175 
SLAY(14)=  0.225 
SLAY(13)=  0.275 
SLAY(12)=  0.325 
SLAY(11)=  0.375 
SLAY(10)=  0.425 
SLAY(9)=  0.497 
SLAY(8)=  0.594 
SLAY(7)=:  0.688 
SLAY(6)=  0.777 
SLAY(5)=  0.856 
SLAY(4)=  0.920 
SLAY(3)=  0.960 
SLAY(2)=  0.981 
SLAY(1)=  0.995 

C  INTERNAL  CONSTANTS 

INSTAB=  0 
ICRILV=  0 
NLEV=  NLAY 
NLEVMU  NLAYMl 
NLEVP1=  NLAYPl 
AKWNMB=  2.5E-05 
LSTAR=  1.0/AKWNMB 
GOCP=  GRAY/ 1005. 

NBASE=  8 

NBASEP1=  NBASE+  1 
K8ASElVil=  NBASE-l 

C  CONSTRAIN  THE  VARIENCE 

III=  MLAT+  1-LATCO 
D0  4  J=  LMLONG 

IF(VAR(J,III).GT.  160000.)  VAR(J,III)=  160000.0 
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4 


CONTINUE 


C  COMPUTE  PRESSURE  AT  EVERY  LAYER 

DO  7  LAY=  1,NLAY 

DO  9  J=  l.MLONG 

PP(J,LAY)=  SLAY(LAY)*PSFC(J) 

9  CONTINUE 

7  CONTINUE 

C  COMPUTE  PRESSURE  AT  EVER  LEVEL 

DO  70  LEV=  1,19 

DO  70  J=  l.MLONG 

PPP(J,LEV)=  SLEV(LEV)*PSFC(J) 

70  CONTINUE 

C  COMPUTE  DENSITY  AT  EVERY  LAYER 

DO  10  LAY=  l.NLAY 

DO  12  J=  l.MLONG 
PRCB=  PP(J,LAY) 

RT=  RGAS*T(J,LAY) 

RO(J,LAY)=  PRCB/RT 

12  CONTINUE 

10  CONTINUE 

C  COMPUTE  DZ  AT  EVERY  LEVEL  FROM  2  TO  18 

DO  14  J=  l.MLONG 
ROILO=  1.0/RO(J.l) 

DO  15  LL=  2,NLEV 
R01UP=  1.0/RO(J,LL) 

R01AVE=  (ROILO+  ROIUP)*0.5 
X=  PP(J,LL-1)-PP(J,LL) 

DZ(J,LL)=  AGRAV*ROIAVE*X 
R01LO=  ROIUP 

IS  CONTINUE 

14  CONTINUE 

C  END  OF  INPUT  AND  ELEMENTARY  COMPUTATIONS 

III=  MLAT^^  1-LATCO 
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DO  17  J=  l.MLONG 


(^***««*«**«,«**««**«*«***««*«****«*«***«*****««**«***«, 

C  SURFACE  AND  BASE  LAYER  STRESS  * 

e****************************************************** 

C  BASE  LAYER  STRESS  IS  DEHNED  IN  TERMS  OF  A  VERTICAL  AVE. 

ROBAR=  0.0 
UBARr:  0.0 
VBAR=  0.0 

DO  200  L=  l.NBASEMl 

C  MASS  WEIGHTED  VERITCAL  AVERAGE  OF  DENSITY,  VELOCITY 

ROBAR=  ROBAR^  RO(J,L)*(PPP(J,L)-PPP(J,L+  1)) 

UBAR=  UBAR4  U(J,L)*(PPP(J,L)-PPP(J,L+  1)) 

VBAR=  VBAR+  V(  J,L)*(PPP(J,L)PPP(J,L+  1)) 

200  CONTINUE 

ROBAR=  ROBAR/(PPP(J,l)PPP(J,NBASE)) 

ROBARs  ROBAR*100.0 

UBAR=  UBAR/(PPP(J,I).PPP(J,NBASE)) 

VBAR=  VBAR/(PPP(J,1).PPP(J,NBASE)) 

C  END  VERTICAL  AVERAGE 

C  VAISALA  FREQUENCY 

DO  201  LEV=  2,NBASE 
LAY=  LEV 

VAI1=  (T(J,LAY)-T(J,LAY1)) 

1  /DZ(J,LEV)+ GOCP 

1F{VAI1.LT.0.0)THEN 
VAI1=  0.0 
ENDIF 

VAI2=  2.0*GRAV/(T(J,LAY) 

1  +  T(J,LAY-1)) 

VSQUA=  VAII*VAI2 
BV(LEV)=  SQRT(VSQUA) 

201  CONTINUE 
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C  VERTICAL  MASS  WEIGHTED  AVERAGE  OF  THE  BRUNT- VAISIALA 
C  FREQ. 

NBAR=  0.0 

DO  202  LEV=  2,NBASE 
LAY=  LEV 

NBAR=  NBAR+  B  V(LEV)*(PP(J,LAY-I)-PP(J,LAY)) 

202  CONTINUE 

NBAR=  NBAR/(PP(J,I)-PP(J,NBASE)) 

IF(NBAR.LE.O.OOOOOOl)  THEN 

PRINT*,’NBAR  IN  GWD  IS  ZERO’ 

ENDIF 

C  DEFINITION  OF  SURFACE  WIND  VECTOR 

UUS=  UBAR 
VVS=  VBAR 

SPEEDS=  SQRT(UUS*UUS+  VVS*VVS) 

IF(SPEEDS.EQ.0.0)THEN 

PRINT*, ’SPEEDS  EQ  ZERO  IN  GWD’ 

ENDIF 

ANG=  ATAN2(VVS,UUS) 

ANGDEG=  180./3.1415926*ANG 

C  STRESS  AT  THE  SURFACE  LEVEL  LEV=  I 

FRSF=  NBAR*SQRT(VAR(J, III))/ SPEEDS 
IF(SPEEDS.EQ.0.0)THEN 
TENSIO(J,I)=  0.0 
ELSEIF(NBAR.NE.0.0)THEN 
GSTAR=  GG(FRSF) 

TENSIO(J,l)=  GSTAR*(ROBAR*SPEEDS*SPEEDS*SPEEDS)/ 

1  (NBAR*LSTAR) 

ELSE 

TENSIO(J,I)=  0.0 
ENDIF 

XTENS(J,I)=  COS(ANG)*TENSIO{J,I) 

YTENS(J,I)=  SIN(ANG)*TENSIO{J,l) 

C  SAVE  SURFACE  VALUES 

DRAGSF(J)=  TENSIO(J,l) 

XDRAG(J)=  XTENS(J.I) 
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YDRAG(J)=  YTENS(J,1) 


ROLO=  RO(J,2) 

UULO=  U(J,2) 

VVLO=  V(J,2) 

TENSIO(J,2)=  TENSICXJ.l) 

XTENS(J,2)=  XTENS(J,1) 

YTENS(J,2)=  YTENS(J,1) 

C  SCALAR  PRODUCT  OF  LOWER  WIND  VECTOR  AND  SURFACE  WIND 

SCALLO=:  UULO*UUS+  VVLO*VVS 

DO  181  LEV=  3,NBASE 
LAY=  LEV 
ROUP=  RO(J,LAY) 

ROAVE=  0.5*(ROLO+  ROUP) 

C  CONVERT  TO  NEWTON/ M2 

ROAVE=  100.0*ROAVE 

C  VELOCITY  COMPONENT  PARALELL  TO  SURFACE  VELOCITY 

UUUP=  U(J,LAY) 

VVUP=  V(J,LAY) 

SCALUP=  UUUP*UUS+  VVUP'VVS 
VELCO=  0.5*(SCALUP+  SCALLO)/ SPEEDS 

C  TAU  DOESN’T  CHANGE  IN  THE  BASE  LAYER  BECAUSE  OF  A 
C  CRITICAL  LEVEL  I.E.  VELCO  <  0.0 

IF(VELCO.LE.O.O)THEN 

TENSIO(J,LEV)=  TENSIO(J,LEV.l) 

GOTO  ISOO 
ENDIF 

C  FROUDE  NUMBER  SQUARED 
YY=  BV(LEV) 

FR02=  YY/(AKWNMB*ROAVE*VELCO*VELCO*VELCO) 

1  •TENSIO(J,LEV.l) 

C  DENOMINATOR  OF  RICHARDSON  NUMBER 
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DELUU=  UUUP-UULO 
DELVV=  VVUP-VVLO 

DELVE2=  (DELLiU*DELUU+  DELVV*DELVV) 


C  RICHARDSON  NUMBER 

IF(DELVE2.NE.0.0)THEN 
DELZ=  DZ(J,LEV) 

VSQUA=  BV(LEV)*BV(LEV) 

RICHSN=  DELZ*DELZ*VSQUA/DELVE2 
ELSE 

RICIISN=  99999.0 
ENDIF 

C  TAU  IN  THE  BASE  LAYER  DOES  NOT  CHANGE  BECAUSE  OF  THE 
C  RICHARDSON  CRITERION 

IF(  RICHSN.LE.0.25)THEN 

TENSIO(J,LEV)=  TENSIO(J, LEV-1) 

GO  TO  1500 
ENDIF 

C  TAU  IN  THE  BASE  LAYER  DOES  CHANGE  IF  THE  LOCAL  FROUDE 
C  EXCEDES  THE  CRITICAL  FROUDE  NUMBER...  THE  SO  CALLED 
C  FROUDE  NUMBER  REDUCTION. 

CRIFRO=  1.0-0.25/ RICHSN 
CRIF2=  CRIFRO'CRIFRO 
IF(LEV.EQ.2)CRIF2=  AMIN  1(0.7, CRIF2) 

IF(  FR02.GT.CRIF2)THEN 

TENSIO(J,LEV)=  CRIF2/FR02*TENSI0(J, LEV-1) 

GO  TO  1500 
ELSE 

TENSIO(J,LEV)=  TENSIO(J, LEV-1) 

GOTO  1500 
ENDIF 

1500  CONTINUE 

XTENS(J,LEV)=  TENSIO(J,LEV)*COS(ANG) 

YTENS(J,LEV)=  TENSIO(  J,LEV)*SIN(ANG) 

ROLO=  ROUP 
UULO=  UUUP 
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VVLO=  VVUP 
SCALLOs  SCALUP 
181  CONTINUE 

C  STRESS  FROM  BASE  LEVEL  TO  TOP  LEVEL 

(;;*«**«*»*««******»«4i*«**«««««***4i**«**«**«**«*********i 


ICRILV=  0 
INSTAB=  0 
ROLO=  RO(J,NBASE) 

UULO=  U(J,NBASE) 

VVLO=  V(J,NBASE) 

C  SCALAR  PRODUCT  OF  LOWER  WIND  VECTOR  AND  SURFACE  WIND 

SCALLO=  UULO*UUS+  VVLO*VVS 

DO  18  LEV=  NBASEP1,NLAY+  1 
LAY=  LEV 

C  THE  STRESS  IS  ALWAYS  INITIALIZED  TO  ZERO 

TENSIO(J,LEV)=  0.0 
IF(ICRILV.EQ.1)THEN 
GO  TO  130 
ENDIF 

1F(LEV.NE.19)THEN 
GO  TO  150 
ENDIF 

TENSIO(J,LEV)=  0.0 
GO  TO  130 

ISO  CONTINUE 

ROUP=  RO(J,LAY) 

ROAVE=  0.5*(ROLO+  ROUP) 

C  CONVERT  TO  NEWTON/ M2 

ROAVE=  100.0*ROAVE 

C  VAIS  ALA  FREQUENCY 

VAI1=  (T(J,LAY)-T(J,LAY-1)) 
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1 


/DZ(J,LEV)+  GOCP 


IF(  VAI1.LT.0.0)THEN 
ICRILV=  1 
TENSIO(J,LEV)=  0.0 
GO  TO  130 
ENDIF 

VAI2=  2.0*GRAV/(T(J,LAY) 

1  +  T(J,LAY-1)) 

VSQUA=  VAI1*VAI2 

C  VAISD  IS  THE  BRUNT- VAISALA  FREQUENCY  N 

VAISD=  SQRT(VSQUA) 

C  VELOCITY  COMPONENT  PARALELL  TO  SURFACE  VELOCITY 

UUUP=  U(J,LAY) 

VVUP=  V(J,LAY) 

C  SCALAR  PRODUCT  OF  UPPER  AND  SURFACE  WIND  VECTOR 

SCALUP=  UUUP*UUS+  VVUP*VVS 
VELCO=  0.5*(SCALUP+  SCALLO)/SPEEDS 
IF(  VELCO.LT.0,0)THEN 
ICRILV=  1 
TENSIO(J,LEV)=  0.0 
GO  TO  130 
ENDIF 

C  FROUDE  NUMBER  SQUARED 
YY=  VAISD 

FR  02=  Y  Y/  ( A  KWN  MB*ROA  VE*  VELCO*VELCO*VELCO) 

C  1  •TENSIO(J, LEV-1) 

C  DENOMINATOR  OF  RICHARDSON  NUMBER 

DELUU=  UUUP-UULO 
DELVV=  VVUP-VVLO 

DELVE2=  (DELUU*DELUU+  DELVV*DELVV) 

C  RICHARDSON  NUMBER 
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IF(  DELVE2.N  E.O.O)TH  EN 
DELZ=  DZ(J,LEV) 

RICHSN=  DELZ*DELZ*VSQUA/DELVE2 
ELSE 

RICHSN=  99999.0 
ICOUNT=  ICOUNT+  1 
ENDIP 

IF(RICHSN.LE.0.2S)THEN 
TENSIO(J,LEV)=  0.0 
ICRILV=  1 
GO  TO  130 
ENDIF 

C  CRITICAL  FROUDE  NUMBER 

CRIFRO=  1.0-0.25/ RICHSN 
CRIF2=  CRIFRO*CRIFRO 

C  END  CRITICAL  FROUDE  NUMBER 

IF(FR02.GE.CRIF2)THEN 
GO  TO  120 
ENDIF 

TENSIO(J,LEV)=  TENSIO(J,LEV-l) 

GO  TO  130 

120  TENSIO(J,LEV)=  CRIF2/FR02*TENSI0(J, LEV-1) 

130  CONTINUE 

XTENS(J,LEV)=  TENSIO(J,LEV)*COS(ANG) 
YTENS(J,LEV)=  TENSIO(J,LEV)*SIN(ANG) 
ROLO=  ROUP 
UULO=  UUUP 
VVLO=  VVUP 
SCALLO=  SCALUP 
18  CONTINUE 

17  CONTINUE 

C  END  STRESS 

DO  20  LAY=  3,NLAY 
DO  21  J=  LMLONG 
LEV=  LAY+  1 
COEF=  GRAV/PSFC(J)*.01 
DSIGMA=  SLEV(LEV)-SLEV(LEV-1) 
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CHUG(J,LAY)=  COKF/DSIGMA*(XTENS(J,LEV)-XTKNS(J, LEV-1)) 
CHVG(J,LAY)=  COEF/DSIGMA*(YTENS(J,LEV)-YTENS(J,LEV-l)) 
21  CONTINUE 

20  CONTINUE 

DSIGMA=  SLEV(3)-SLEV(1) 

DO  444  J=  l.MLONG 

COEF=  GRAV/PSFC(J)*.01 

CHUG(J,1)=  COEF/DSIGMA*(XTENS(J,3)-XTENS(J,l)) 
CHVG(J,1)=  COEF/DSIGMA*(YTENS(J,3)-YTENS(J,l)) 
CHUG(J,2)=  CHUG{J,1) 

CHVG(J,2)=  CHVG(J,1) 

444  CONTINUE 

III=  MLAT+  1-LATCO 

C  ZNORM=  #  OF  DAYS  •  i  OF  STEPS  PER  DAY 

ZNORM=  10.0*96.0 

DO  333  J=  1.96 
DO  333  K=  1,18 

UTEN(J.III,K)=  UTEN(J,III,K)+  CHUG(J,K)/ZNORM 
VTEN(J,III,K)=  VTEN(J,III,K)+  CHyG(J,K)/ZNORM 

333  CONTINUE 

DO  334  J=  1,96 

TAUSX(J.III)=  TAUSX(J,III)+  XDRAG(J)/ZNORM 
TAUSY(J,III)=  TAUSY(J,III)+  YDRAG(J)/ZNORM 

334  CONTINUE 

RETURN 

END 

FUNCTION  GG(FR) 

REAL  FR 
G=  1.0 
A=  1.0 

GG=  (FR*FR)/(FR*FR+  A*A) 

RETURN 

END 
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